package cs.ubbcluj.ro.linearsystem.implementations;

import cs.ubbcluj.ro.datatypes.interfaces.Matrix;
import cs.ubbcluj.ro.datatypes.interfaces.Vector;
import cs.ubbcluj.ro.linearsystem.interfaces.LinearSystem;

public class LinearSystemLU implements LinearSystem {
	
	Matrix A;
	Vector b;
	
	public LinearSystemLU() {
	}

	public LinearSystemLU(Matrix A, Vector b) {		
		this.A = A;
		this.b = b;
	}
	
	/**
	 * calculates the L matrix given by the Cholesky algorithm
	 * @param m
	 * @return
	 * @throws InterruptedException 
	 */
	Matrix calculateLMatrix(Matrix m) throws InterruptedException {
		// create Helper Matrix to store the result
		Matrix resultMatrix = m.createHelperMatrix(m.getNrOfRows(), m.getNrOfColumns());
		
		
//		System.out.println("\n-----------------------calculateLMatrix---------------------");
//		
//		System.out.println("\nm = ");
//		m.print();
//		
//		System.out.println("\nresultMatrix = ");
//		resultMatrix.print();
		
		// initialize sum variable
        double sum;
        
        // row indices
        for (int k = 0; k < m.getNrOfRows(); k++) {	
//        	System.out.println("k = " + k);
        	// set diagonal elements
            sum = 0;           
            for (int j = 0; j < k; j++) {
//            	System.out.println("j = " + j);
//            	System.out.println("resultMatrix[k,j] = " + resultMatrix.get(k, j));
                sum += Math.pow(resultMatrix.get(k, j),2);
//                System.out.println("sum = " + sum);
            }
            resultMatrix.set(k, k, Math.sqrt(m.get(k,k)-sum));
            
//            System.out.println("Math.sqrt(m.get(k,k)-sum) = " + m.get(k,k));
//            System.out.println("resultMatrix[k,k] = " + resultMatrix.get(k, k));
            
            // column indices
            for (int i = k+1; i < m.getNrOfRows(); i++) {
//            	System.out.println("i = " + i);
                sum = 0;
                for (int p = 0; p < k; p++) {
//                	System.out.println("p = " + p);
                    sum += resultMatrix.get(i,p) * resultMatrix.get(k,p);
 //                   System.out.println("sum = " + sum);
                }
                resultMatrix.set(i,k,((m.get(i,k)-sum)/resultMatrix.get(k,k)));
            }            
        }
        
//        System.out.println("-----------------------calculateLMatrix---------------------");
        return resultMatrix;
    }
	
	/**
     * forward calculation
     * calculates the vector z given by the formula: Lz = b, where L is the lower matrix of A
	 * @param m
	 * @param v
	 * @return
	 * @throws InterruptedException
	 */
    Vector forwardSolveCycle(Matrix m, Vector v) throws InterruptedException {
    	Vector resultVector = v.createHelperVector(v.getNrOfRows());
        double sum;
        for (int k = 0; k < v.getNrOfRows(); k++) {
            sum = 0;            
            for (int p = 0; p < k; p++)
                sum += resultVector.get(p)*m.get(k,p);
            resultVector.set(k,(v.get(k)-sum)/m.get(k,k));
        }
        return resultVector;
    }
    
    /**
     * backward calculation
     * calculates the vector z given by the formula: L^T * x = z, where L is the lower matrix of A
     * @param m
     * @param v
     * @return
     * @throws InterruptedException
     */
    Vector backwardSolveCycle(Matrix m,Vector v) throws InterruptedException {
    	Vector resultVector = v.createHelperVector(v.getNrOfRows());    	
        double sum;
        
        for(int k = m.getNrOfRows()-1; k>=0; k--) {
            sum = 0;            
            for(int p = k+1; p < m.getNrOfRows();p++)
                sum += m.get(k,p)*resultVector.get(p);
            resultVector.set(k,(v.get(k)-sum)/m.get(k,k));
         }
        return resultVector;
    }

	/**
	 * Solves the linear system using Cholesky decomposition.
	 * Formula: Ax = b
	 */ 
	public Vector solveSystem() throws InterruptedException {		
		
		// decomposition
		Matrix lMatrix = calculateLMatrix(A);
		
//		System.out.println("\nlmatrix");
//		lMatrix.print();
		
		// forward solving 
		Vector z = forwardSolveCycle(lMatrix, b);
		
//		System.out.println("\nz");
//		z.print();
		
		// backward calculation
		Vector resultVector = backwardSolveCycle(lMatrix.trans(), z);
		
//		System.out.println("\nresutVector");
//		resultVector.print();
		
		return resultVector;
	}
	
	public Vector solveSystem(Matrix A, Vector b) throws InterruptedException {
		// set variables
		this.A = A;
		this.b = b;
		
		// solve system
		return solveSystem();
	}
}
